function [delta,invW,fx1,fz,fyp] = nextStepFixedPointSolverBoucekkine(Z,x_t,y_tN,N,...
    f,invW,fx1,fz,fyp,updateJ,setupEPer)
% This function solves the linear system for the Newton-Raphson fixed-point
% algoirithm using the algorithm outlined in Boucekkine (1995, JEDC)
% Notation:
% The matrices R and Q, are structured as follows for a given time period:
% R(:,:,k) = [R_ny; R_nx1]
% Q(:,:,k) = [Q_ny; Q_nx1]

% Unfold setupEPer
nx        = setupEPer.nx;
mx        = setupEPer.mx;
myx       = setupEPer.myx;
ny        = setupEPer.ny;
params    = setupEPer.params;
h0        = setupEPer.h0;
hx        = setupEPer.hx;
g0        = setupEPer.g0;

[x,y]    = VarStack(Z,nx,mx,myx,ny,hx,x_t,y_tN,N);
R        = zeros(ny+mx,1,N);
Q        = zeros(ny+mx,ny,N);

if updateJ == 1
    getDeriv = 1;
    invW     = zeros(ny+mx,ny+mx,N);
    fx1      = zeros(ny+mx,mx,N);
    fyp      = zeros(ny+mx,ny,N);
    fz       = zeros(ny+mx,myx,N);
    % First period
    k = 1;
    [~,fx1(:,:,k),fx1p,fy,fyp(:,:,k),fz(:,:,k)] = getModelDeriv_EP(nx,mx,myx,ny,params,N,h0,g0,x,y,k,getDeriv);        
    invW(:,:,k) = [fy, fx1p]\eye(ny+mx);
    Q(:,:,k)    =  invW(:,:,k)*fyp(:,:,k);
    R(:,:,k)    = -invW(:,:,k)*(fx1(:,:,k)*x_t(1:mx,1)+f(:,k));
    % Remaining periods
    for k=2:N
        [~,fx1(:,:,k),fx1p,fy,fyp(:,:,k),fz(:,:,k)] = getModelDeriv_EP(nx,mx,myx,ny,params,N,h0,g0,x,y,k,getDeriv);        
        invW(:,:,k) = [-fz(:,:,k)*Q(1:myx,:,k-1)-fx1(:,:,k)*Q(ny+1:ny+mx,:,k-1)+fy, fx1p]\eye(ny+mx);
        Q(:,:,k)    =  invW(:,:,k)*fyp(:,:,k);
        R(:,:,k)    = -invW(:,:,k)*(fz(:,:,k)*R(1:myx,:,k-1)+fx1(:,:,k)*R(ny+1:ny+mx,:,k-1)+f(:,k));
    end
else
    % First period 
    k = 1;
    Q(:,:,k) =  invW(:,:,k)*fyp(:,:,k);
    R(:,:,k) = -invW(:,:,k)*(fx1(:,:,k)*x_t(1:mx,1)+f(:,k));
    % Remaining period
    for k=2:N
        Q(:,:,k) =  invW(:,:,k)*fyp(:,:,k);
        R(:,:,k) = -invW(:,:,k)*(fz(:,:,k)*R(1:myx,:,k-1)+fx1(:,:,k)*R(ny+1:ny+mx,:,k-1)+f(:,k));
    end
end
ySolv  = zeros(ny,N);   %containing y_t, y_t+1,...,y_t+N-1
x1Solv = zeros(mx,N);  %containing x1_t+1,x1_t+2,...,x1_t+N
% Last period
k = N;
ySolv(:,k)  = R(1:ny,:,k)-Q(1:ny,:,k)*y_tN;
x1Solv(:,k) = R(ny+1:ny+mx,:,k)-Q(ny+1:ny+mx,:,k)*y_tN;
% Remaining periods
for k=N-1:-1:1
    ySolv(:,k)  = R(1:ny,:,k)       - Q(1:ny,:,k)*ySolv(:,k+1);
    x1Solv(:,k) = R(ny+1:ny+mx,:,k)- Q(ny+1:ny+mx,:,k)*ySolv(:,k+1);
end
% The solution to the linear system in the Newton-Raphson routine
delta = reshape([ySolv;x1Solv],(ny+mx)*N,1);

end

